↜ Back to index Introduction to Numerical Analysis 2

Solutions for Lecture 6

Exercise 3

implicit none
integer n, i, j, k
real, dimension(:,:), allocatable :: a
real, dimension(:), allocatable :: b, x
real c, s

print *, 'Enter n, a, b'
read *, n

allocate(a(n, n))
allocate(b(n))
allocate(x(n))

do i = 1,n
    read *, a(i, :)
end do

read *, b

! elimination step
do k = 1,n-1
    ! add k-th row to the remaining ones
    do i = k+1,n
        if (a(k, k) == 0) then
            print *, 'Error: a(', k, ',', k, ') is zero'
            stop
        end if
        c =  a(i, k) / a(k, k)
        do j = k,n
            a(i, j) = a(i, j) - a(k, j) * c
        end do 
        ! fancy way:
        ! a(i, :) = a(i, :) - a(k, :) * c
        b(i) = b(i) - b(k) * c
    end do
end do

do i = n,1,-1
    s = 0.
    do j = i + 1, n
        s = s + x(j) * a(i, j)
    end do 
    x(i) = (b(i) - s) / a(i, i)
end do

print *, 'x = '
do i = 1,n
    print *, x(i)
end do

! print *, 'a = '
! do i = 1,n 
!     print *, a(i, :)
! end do
!
! print *, 'b = '
! do i = 1,n
!     print *, b(i)
! end do

end